using System;
using System.Reflection;
using NBody;
using NBody.Analysis;
using NBody.Cosmology;
using OptionParser;
using DataFile;

[assembly: AssemblyTitle ("NBodyBar")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Calculate the bar strength of an N-body system.")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class BarOptions : InputOutputGetOptions
{
    [Option ("-c", "Adjust for centre of mass")]
    public bool CentreOfMass;
    
    [Option ("-l", "Exclude particles lower than this radius.")]
    public double LowerLimit;
    
    [Option ("-u", "Exclude particles greater than this radius.")]
    public double UpperLimit;
    
    [Option ("-v", "Exclude particles outside virial radius.")]
    public bool ExcludeVirial;
    
    [Option ("-cosmo", "Cosmology for virial radius calculation.")]
    public string CosmoString;
    
    [Option ("-b", "Calculate bar strength at \"b\" intervals in radius.")]
    public int NumBins;
    
    [Option ("-e", "Calculate bar strength from all particles enclosed in radius,\n              rather than in radial bins.")]
    public bool Enclose;
    
    [Option ("-sl", "Semi-log plot.")]
    public bool SemiLog;
    
    public BarOptions()
    {
        LowerLimit = 0.0;
        UpperLimit = Double.PositiveInfinity;
        CosmoString = "SCDM";
        NumBins = 0;
    }
}

public class NBodyBar
{
    public static void Main(string[] args)
    {
        // Process command-line options
        BarOptions opts = new BarOptions();
        opts.ProcessArgs(args, "bar");
                        
        // Read in a N-body system from the given file
        NBodySystem s = NBodySystem.Read(opts.InFile);
        
        // Right now, the user can specifiy to adjust for c/m or not, 
        // but not the tolerance parameter.  I think this is fine.
        if (opts.CentreOfMass)
            s.AdjustForCentreOfMass(1e-3);
        
        // if the user wishes to exclude particles outside the virial
        // radius, we need a cosmology to calculate the background density
        // at s.Time.  The user can input the cosmology as a string, but
        // it better corresponde to a specific cosmology defined in the 
        // Cosmology class.
        if (opts.ExcludeVirial)
        {
            Model cosmo = Model.FromString(opts.CosmoString);
            double z = cosmo.Redshift(s.Time);
            opts.UpperLimit = VirialRadius.CalculateR200(s, cosmo.Density(z));
        }
        
        // If the user didn't specify the number of bins, just compute
        // the bar strength over the whole system
        if (opts.NumBins == 0)
        {
            double bar_strength = Bar.CalculateStrength(s, opts.LowerLimit, opts.UpperLimit);
        
            Console.Error.Write("Time: ");
            Console.WriteLine(s.Time);
            Console.Error.Write("Bar strength: ");
            Console.WriteLine(bar_strength);
        }
        // But if he did, we need to either:  bin the system and calculate
        // the bar strength in each bin, or calculate the bar strength
        // from all the particles up to the bin radius.
        else
        {
            Table data = new Table(opts.NumBins, 2, "Radius, Bar Strength");
            if (Double.IsPositiveInfinity(opts.UpperLimit))
                opts.UpperLimit = s.Rmax;
            
            for (int i = 0; i < opts.NumBins; i++)
            {
                data[i, 0] = opts.LowerLimit + (double)(i+1) / (double)opts.NumBins * (opts.UpperLimit - opts.LowerLimit);
                if (opts.Enclose)
                    data[i, 1] = Bar.CalculateStrength(s, opts.LowerLimit, data[i, 0]);
                else if (i == 0)
                    data[i, 1] = Bar.CalculateStrength(s, opts.LowerLimit, data[i, 0]);
                else
                    data[i, 1] = Bar.CalculateStrength(s, data[i-1, 0], data[i, 0]);
            }
            
            // done; print it out.  Can also print with log/log, but ....
            if (opts.SemiLog)
                data.Print(opts.OutFile, new Transform[] { Math.Log10, Table.IdentityTransform });
            else
                data.Print(opts.OutFile);
        }                
    }
}
